Quantum Monte Carlo Study of the Disordered Attractive Hubbard Model 
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We investigate the disorder-driven superconductor to insulator quantum phase transition (SIT) in 
an interacting fermion model using determinantal quantum Monte Carlo (QMC) methods. The 
disordered superconductor is modeled by an attractive Hubbard model with site disorder chosen 
randomly from a uniform distribution. The superconducting state which exists for small disorder is 
shown to evolve into an insulating phase beyond a critical disorder. The transition is tracked by the 
vanishing of (a) the superfluid stiffness, and (b) the charge stiffness or the delta function peak in 
the optical conductivity at zero frequency. We also show the behavior of the charge, spin, pair, and 
current correlations in the presence of disorder. Results for the temperature dependence of the dc 
conductivity, obtained by an approximate analytic continuation technique, are also presented both 
in the metallic phase above T c and the insulating phase. We discuss some of the complications in 
extracting the resistance at the transition point. 
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I. INTRODUCTION 



In a wide variety of two dimensional disordered sys- 
tems jjj, from granular and homogeneously disordered 
Bi, Pb and Sn films to Im^CX, § and MoGe 

films ||, high temperature superconducting films [0J|] 
and Josephson junction arrays a transition from a su- 
perconductor to an insulator (SIT) can be driven by ad- 
justing some tuning parameter such as the film thickness, 
the O concentration, the magnetic field strength, or the 
charging energy. The experimental signature of the tran- 
sition is that the behavior of the sheet resistance Rn (T) 
as a function of temperature T is different in the two 
phases. At low disorder or magnetic field, the system is 
superconducting for T < T c . The transition temperature 
T c decreases with increasing disorder or magnetic field 
and above T c the system is metallic with dRa/dT > 0. 
Beyond a critical disorder or magnetic field, on the other 
hand, the system becomes insulating with dR n /dT < 0. 

Motivated by these experiments, one of the important 
open theoretical questions is to study particular micro- 
scopic models to see whether or not they show a SIT as 
a function of some tuning parameter such as the degree 
of disorder and, if so, characterize the transition. 

Anderson [Tc| ] proposed that the superconducting tran- 
sition temperature T c and the thermodynamic properties 
should be unaffected by disorder since Cooper pairs can 
be formed by pairing the time-reversed exact eigenstates 
of the noninteracting disordered problem. This is only 
valid for small disorder in the regime hp I 3> 1, where 
kp is the Fermi momentum and t is the elastic mean 
free path. Ma and Lee |ll| developed a mean field the- 
ory in which they assumed that the order parameter was 
uniform throughout the system. As a consequence, the 
superfluid density remained large even for fairly high dis- 



order and was found to persist essentially all the way to 
the site localized limit. 

One might therefore ask whether a disorder-driven SIT 
can occur at all. It is important to note that both the An- 
derson and Ma-Lee arguments make specific assumptions 
concerning the effect of randomness, and hence may not 
be compelling in all cases. In order to understand why 
a SIT might be possible, consider the two generic mech- 
anisms for the destruction of superconductivity. First, 
the magnitude of the pairing gap can be driven to zero. 
Second, phase coherence between the pairs in different 
parts of the sample may be lost. Clearly there is an in- 
terplay between fluctuations in the pair amplitude and 
phase. For example, the phase can change at a smaller 
energy cost in regions where the amplitude is lower |L2| . 
It is possible that the pair amplitude is driven to zero at 
the same point where phase coherence is lost, but it is 
also possible that the two phenomena occur separately. 

Fisher and collaborators Jl3| were the first to describe a 
scenario in which phase fluctuations caused a SIT while 
the pair amplitude remained finite. They conjectured 
that the SIT might be in the same universality class 
as the superfluid-insulator transition for bosons. They 
argued that since near the transition the size of the 
Cooper pair is much smaller than the diverging correla- 
tion length, it is possible to describe it as a bose field. Of 
course, the charge carriers of the experimental systems 
are fermionic in nature, so it is useful to study Hamilto- 
nians which do not begin immediately with bosonic de- 
grees of freedom. Perturbative methods to study the SIT 
in fermionic models have not been successful in describ- 
ing the transition region fl^ , [f5| , which is not surprising 
since the transition occurs in a region of high disorder in 
an interacting system. 

While this approach has led to a number of very inter- 
esting results, especially for the value of the conductiv- 
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ity at the transition [|l6|-|l9|] , it is important to test the 
validity of the phase-only models by developing meth- 
ods which also treat amplitude fluctuations. In order to 
better describe the behavior of a superconductor at high 
disorder, Ghosal and collaborators pQ] have included the 
fluctuations of the superconducting order parameter by 
solving the "Bogoliubov-de Gennes" mean field equations 
self-consistently. They have found that the probability 
distribution of the local pairing amplitude develops a 
broad distribution with significant weight near zero with 
increasing disorder. Surprisingly, the density of states 
continues to show a finite spectral gap, as also seen by 
Quantum Monte Carlo (QMC) and maximum entropy 
techniques [|T), shown to arise from the break up of 
the system into superconducting islands separated by re- 
gions with very small pairing amplitude. These disorder- 
induced fluctuations in the order parameter amplitude 
have a marked effect in suppressing the superfluid den- 
sity at higher disorder but by themselves are not suf- 
ficient to drive the system non-superconducting. It is 
necessary to include phase fluctuations distributed inho- 
mogeneously riding on top of the highly inhomogeneous 
amplitude fluctuations to get a SIT. 

In this paper we describe the first QMC study of a 
fermion model of superconductivity (the attractive Hub- 
bard Hamiltonian with random site energies) which gives 
a SIT at a critical disorder strength p2[ . The attrac- 
tive Hubbard Hamiltonian which we study is a simple 
model of a disordered SC that allows us to explore the 
qualitative issues arising from the interplay of supercon- 
ductivity and localization. While such a model does not 
address questions concerning the microscopic origin of 
the pairing, since the attraction is put in a priori, one 
can nevertheless examine questions such as the competi- 
tion between superconductivity and charge density wave 
formation p3| , the behavior of superconducting corre- 
lations above the superconducting transition tempera- 
ture p4|-|2"8|] , and the interpolation between weak cou- 
pling BCS and strong coupling bosonic regimes of pair 
formation 29 1. 



II. ORGANIZATION OF THE PAPER 

This paper is organized as follows: In section III we in- 
troduce the attractive Hubbard model and briefly review 
the physics of the clean attractive Hubbard model. In 
Section IV we describe the QMC simulation technique. 
In section V we first discuss our results for the chem- 
ical potential, in order to demonstrate that we are in 
the degenerate Fermi regime of the model. We then 
describe the effect of disorder on the local and longer 
range density-density and pairing correlations. The pair- 
ing correlations are found to be much more robust com- 
pared to the density correlations away from half filling. 
We also show the behavior of the superconducting order 
parameter which decreases rapidly with increasing disor- 



der and vanishes beyond a critical disorder. In section 
VI we present a detailed discussion of the longitudinal 
and transverse current-current correlation functions. The 
longitudinal response obeys the f-sum rule and equals the 
absolute value of half the lattice kinetic energy K x which 
we verify in our simulations. The transverse response on 
the other hand, deviates from K x and this deviation is 
a measure of the superfluid stiffness of the system. We 
present results showing the suppression of the superfluid 
stiffness with disorder and its ultimate destruction be- 
yond a critical disorder. In section VII we discuss the 
behavior of the frequency dependent current-current cor- 
relation function and the extraction of the charge stiffness 
or the strength of the delta-function peak in the optical 
conductivity. Our results show that in the superconduct- 
ing phase, the superfluid stiffness and the charge stiffness 
are roughly equal in magnitude for all disorder strengths. 
In section VIII we discuss an approximate method to ex- 
tract the temperature dependence of the dc resistivity 
and show its behavior in the metallic phase above T c for 
low disorder as well as in the insulating phase for higher 
disorder. The resistivity at the transition is extracted by 
two methods- (i) At the critical disorder, the charge stiff- 
ness vanishes with frequency with a slope proportional to 
the resistivity; and (ii) from the crossing of the resistivity 
vs disorder curves at various temperatures. We also dis- 
cuss the complications of obtaining the resistivity near 
a quantum critical point. We present our conclusions in 
section IX and end with some of the outstanding ques- 
tions in the area of SIT in section X. In previous papers 
p2| , p0 31 1 we have presented a short discussion of some 
of these issues. The purpose of the present manuscript 
is to provide the details behind that work, as well as to 
present a number of new results including a more com- 
plete discussion of both the physics and the numerics. 



III. MODEL 



The Hamiltonian we study is defined by 



H = I>lcja + c} CT C iCT ) 



■5>"Vi)n iCT -|t/|E^T-|)Kl-|) ' (1) 



Here the lattice sum (ij) is over nearest neighbor sites on 
a two dimensional square lattice, c- lr7 is a fermion destruc- 
tion operator at site i with spin er, n- la = ^ G c\ a , and the 
chemical potential fi fixes the average density (n). The 
site energies V\ are independent random variables with a 
uniform distribution over [— V/2,V/2]. The interaction 
has been written in particle-hole symmetric form so that 
/j = corresponds to (n) — 1 at all U and T when V = 0. 
We set t = 1 and measure all energies in units of t. 

In real materials, disorder plays a complicated role 
in the Hamiltonian, both affecting the screening of the 
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electron-electron interaction as well as the phonons and 
hence the electron-phonon interaction. Our Hamiltonian 
does not include these effects. 

Some of the physics of the clean attractive Hubbard 
model may be summarized as follows p2| , ^3[ : At half- 
filling (/i = 0), the model has no long range correla- 
tions at any finite temperature, and at T = is in a 
state with combined charge density wave (cdw) and su- 
perconducting order [ p4[ . When fj, =/= the system has 
a finite temperature Kosterlitz-Thouless transition to a 
state with superconducting order. The transition tem- 
perature T c depends strongly on the filling near (n) = 1 . 
T c shows a non-monotonic dependence on coupling |2^] 
similar to the repulsive Hubbard model where the Neel 
temperature first increases with U but then goes down 
as Tjy oc J = At 2 /U at strong coupling. Numerical es- 
timates | f}5|]3"6l of T c are still a matter of considerable 
debate and at (n) = 0.875 vary from 0.3t to 0.03i. 



IV. QUANTUM MONTE CARLO SIMULATION 



sum over all configurations of the HS fields on the dis- 
cretized space-time lattice. The partition function in the 
grand canonical ensemble is 



Z = Trcxp 



exp 



{■5} 



C h h {S}(T, v)Cjc 



(•5) 



Here h{ $} (t, o~) is a one-body Hamiltonian for the motion 
of an electron in a given configuration of the H-S fields. 
Note in Eq. || both the up and down electrons couple to 
the HS field with the same sign. 

Now the resulting trace over quadratic forms in the 
fermion operators in Eq. |5| is performed and gives 



5^detM T ({5})detM;({5}) 

{S} 



(6) 



with 



Our simulation uses the standard "determinant" QMC 
algorithm along with its various refinements 

@|o|. The partition function Z = Tr[e _/3ff ] is writ- 
ten as a path integral by discretizing the imaginary time 
dimension = 1 /T into N T time slices as 



-/3H 



AtH\ n t 



(e -ArH le -Ar Hu) N, 



(2) 



where [3 — N t At. In Eq. 0, Hi is the sum of the two 
single particle terms in Eq.^ and Hjj is the interaction 
term. A systematic Trotter error is introduced in Eq. |^ 
because of the non-commutativity of the operators Hi 
and Hjj. This Trotter error, however, can be dealt with, 
either by making At sufficiently small so that errors in 
observables are of the same order as statistical fluctua- 
tions from the sampling, or, if greater accuracy is needed, 
by extrapolating to At = 0. The exponential of the inter- 
action term is decoupled using a Hubbard-Stratonovich 
(HS) transformation by introducing a discrete field 0] 
Si T = ±1 at each point in the space-time lattice, 

exp[+AT | U | K T - l/2)(nu - 1/2)] = 



2 exp 
where 



{ AT } Ul } E exp[ATAS lT (n lT +n a -l)] (3) 
L J S w =±l 



cosh(ATA) = cx P (At | U | /2), 



(4) 



is satisfied by real A. Thus the original functional integral 
over Grassman variables, which involved traces contain- 
ing quartic operators is reduced to a quadratic problem 
in the fermion operators but at the cost of performing a 



M a {{S}) 



-irJi( S )(T,j) 



(7) 



Thus the interacting problem is equivalent to solving a 
non-interacting problem for a given HS field configura- 
tion {Si T } and then summing over all possible configu- 
rations. The sum over the HS fields on the space-time 
lattice is efficiently done using Monte Carlo techniques 
which generate the configurations, treating the product 
of the determinants as a probability. Note that in general 
for a fermion problem, since the sign of the determinants 
may be negative, the product is not necessarily non- neg- 
ative and it cannot be treated as a probability. This is 
the origin of the 'sign-problem' for typical fermion prob- 
lems. However, for the Hamiltonian in Eq. [j], since it is 
possible to couple the HS field to the charge and 
satisfy Eq. || with real A, the two determinants in Eq. || 
are identical, and hence the integrand is non-negative- 
thus there is no sign problem |41| in attractive Hubbard 
model simulations at any filling. 

In the determinant QMC approach, finite temperature 
expectation values of combinations of fermion operators 
with arbitrary space and imaginary time arguments can 
be easily evaluated. More precisely, if all the operators 
are at the same imaginary time, the observables can be 
expressed in terms of matrix elements of the inverse of the 
matrices whose determinants give the Boltzmann weight. 
These matrix elements are needed to update the HS field, 
and are therefore available "free of charge" for the mea- 
surements. If the operators whose expectation values are 
to be measured have different imaginary time arguments, 
some extra calculations are involved to obtain the non- 
equal time Green's functions. However this can be done 
in a straightforward manner Hullo)- 
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V. EQUAL TIME CORRELATIONS 



B. Density-density correlations 



A. Chemical Potential 

The location of the chemical potential relative to the 
bottom of the band gives information about the degen- 
eracy of the system. In the simulations presented in this 
paper the filling is chosen to be (n) = 0.875 close to the 
poi nt where T c is expected to be maximal for U = — At 
p5[ . For a given value of the parameters-interaction 
strength U, disorder strength V and temperature T-the 
chemical potential is tuned so that upon disorder av- 
eraging the density (n) ~ 0.875. We comment that an 
alternative approach is to tune the chemical potential for 
each disorder realization separately so that each has the 
same desired filling. This is likely to result in reduced 
fluctuations f42|| , but is considerably more time consum- 
ing numerically. Some such approach, however ; appears 
essential for analytic continuation calculations |21| . 
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FIG. 1. The chemical potential /i shows a roughly lin- 
ear decrease with disorder V. Since (j,(T,\U\,V) + At+ 
{n)\U\/2 ~ 3.5i S> T the system is in a highly degenerate 
regime and away from the preformed bosonic regime. 



In Fig. |^ we show the double occupancy (njfiiij.) which 
is found to increase from 0.32 at V = to 0.38 at V — 5. 
This increase is a consequence of the fact that in the 
attractive model, random site energies and interactions 
both act to promote double occupancy, in contrast to the 
repulsive model where they compete. 




FIG. 2. The increase in double occupancy (nifnij.) as a 
function of increasing disorder V/t. 

In Figs. ||(a) and (b) we also show the spatial variation 
of the density- density correlation function 



Ca,<r'(l) = (n iCT n i+ i i(T /} - (n iCT )(n i+ i lCT ») 



(8) 



At half filling, C(l) is rapidly suppressed by disorder p3| ; 
via finite size scaling it is seen that even as little disorder 
as V = Q.25t is capable of destroying the charge density 
wave ordering and in an 8 x 8 system C(l) is definitely 
suppressed by V — It. Away from half filling even for the 
clean system C(l) is small and thereafter disorder does 
not have any further effect. 



The dependence of fi on V is roughly linear and is 
shown in Fig. [I] for U = —At. Since fj,, measured from the 
bottom of the band and taking into account the Hartree 
shift, is larger than the temperature, /x(T, \ U\, V) + At + 
(n)\U\/2 > T the system is degenerate and far from the 
regime where there are preformed bosons. Note, we have 
assumed that the bottom of the band is at — At, which is 
the case in the clean system but should be renormalized 
by the random potential in the disordered system. 



C. Pair Correlations 

An important characteristic of the superconducting 
state is that the equal time s-wave pair correlation func- 
tion P s defined by, 

p s {\) = ( a-aU >, 

At = c t TC t ;) ( 9 ) 

has a finite value at large separations P s (l = 
(i/2,L/2)) = Aq P , where Aop is the "order param- 
eter" on a lattice of finite size L. 
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FIG. 3. The density correlation function C CTjCr /(l) from 
Eq. | for 1 along [10] and [11] directions for (a) (a, <j') = (f, |) 
and (b) (a, a') = (T,T), showing rapid suppression with in- 
creasing disorder. 

In Fig. [| we show the behavior of P s at a temper- 
ature T = O.lt for varying degrees of disorder. This 
temperature is sufficiently low that for the clean system 
the correlation length has exceeded the linear lattice size 
and the system is effectively in the ground state. For 
the clean system, or weak disorder, the correlation func- 
tion approaches a constant at large distances, implying 
a SC state with long range order. For strong disorder, 
the correlation function vanishes at large distances indi- 
cating the absence of an order parameter. It is evident 
by comparison with Fig. || that pairing correlations are 
much more robust than density-density correlations for 
the same degree of disorder, as in the half filled case P3[. 
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FIG. 4. The pair correlation function defined in Eq. [] 
is shown as a function of 1, the relative separation of the 
two sites along [10] and [11] directions for varying disorder 
strengths V =0, 1.0. 2.0, 3.0, 3.5, 4.0, and 5.0. The value at 
1 = is given by Eq. [To] but is not shown as it is off-scale. Note 
the relative robustness of the pairing correlations compared to 
the density correlations in Fig. |§] in the presence of disorder. 
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FIG. 5. Suppression of the superconducting "order param- 
eter" Aop on an 8x8 lattice with increasing disorder. While 
Aop does not vanish at large V due to finite size effects, 
a scaling analysis of the pair structure factor indicates that 
in the thermodynamic limit Aop vanishes around a critical 
disorder V c ~ 3.5t. 

Fig. ^ shows the order parameter Aop as a function 
of disorder which is strongly suppressed by disorder and 
vanishes beyond a critical disorder strength V c ~ 3.5i. 

The value of the pairing correlation function at zero 
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separation is related to the occupancy and double occu- 
pancy, 

P s (0) = (A,Aj) 

= l - (n) + (mimi) • (io) 

Whereas P s (l) is reduced by disorder for 1 nonzero, P s (0) 
is increased, since the density (n) is fixed and the double 
occupancy rate (n^n^) is increased (Fig. 2). 



superconductivity. The long range pairing order in the 
ground state is suppressed to zero for disorder V ~ At, 
when U — —At. Off half-filling, the charge correlations 
are small and little affected by randomness, though disor- 
der does cause an enhancement of the double occupancy 
rate. However, considerably more information can be ob- 
tained by looking also at various imaginary time depen- 
dent quantities such as the current-current correlation 
function. 



00 
6 



- d 



o 



n i i | i i i | i i i | r 

U = -4 V=1 (°) 







oq 
d 



- d 



o 







aT=0.50 
□ T=0.17 
oT=0.10 



J I I I I L 



_L 



_l I I I L 



1 2 



3 



i — i — i — | — i — i — i — | — i — i — i — | — r 

U = -4 V=4 ( b ) 



aT=0.50 
□ T=0.17 
oT=0.10 

J I I I I L 



1 



2 



3 



FIG. 6. The longitudinal current-current correlation func- 
tion A xx (q x ) defined in Eq. |l2| as a function of q x at T = 0.5t 
(open triangles), 0.17i (open squares), and O.lt (open circles). 
The corresponding filled points at q x = are the magnitude 
of the kinetic energy K x along x at those temperatures. In 
(a) V = It and in (b) V = At. In all cases A L = A xx (q x -» 0) 
approaches K x as required by gauge invariance. 

The equal time pair and density correlations already 
give considerable insight into the effect of disorder on 



VI. CURRENT-CURRENT CORRELATION 
FUNCTION 

As known for some time j43|, and also described re- 
cently in the context of quantum simulations f44|| , various 
limits of the current-current correlation function give in- 
formation about the charge and superfluid stiffness, and 
gauge invariance, and in principle can be used to dis- 
tinguish insulators, metals, and superconductors. The 
current-current correlation function A xx (1,t) is defined 
by 



A M (l,r) = (j x (l,T)j x (0,0)) 



j»(lT) = e HT 



HT (11) 



Upon Fourier transforming in space and imaginary time 
we get A xx (q,uj n ) = J2i /<f dre^e' 1 ^ A xx (l, r), where 
u) n = 2nn/p. 



A. Longitudinal response 

The longitudinal part of A xx defined in Eq. [ll] must 
satisfy the f-sum rule, 

A L = lim 9x ^ A xx (q x ,q y = 0, ut n = 0) 

A L = K x (12) 

as a consequence of gauge invariance . Here K x = 

(t Y,A4+x,v c h<r + c {a c \+x.a)) is the magnitude of the ki- 
netic energy in the x direction. 

Fig. ^| shows A xx {q x ) as a function of q x for different 
temperatures at weak disorder V = It (in (a)) and at 
strong disorder V = At (in (b)). In both cases one finds 
that A L = A xx (q x — > 0) = K x at all T, verifying the 
gauge invariance condition and providing a non-trivial 
check of our numerics. 



B. Transverse response: Superfluid Stiffness 

The transverse response is given by 

A T = linig y ^o A xx (q x = 0, q y ,uj n = 0) • (13) 
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FIG. 7. The transverse current-current correlation func- 
tion A X x(q y ) defined in Eq. [l3] as a function of q y at T = 0.5t 
(open triangles), 0.17t (open squares), and O.lt (open circles). 
The corresponding filled points at q y — are the magnitude 
of the kinetic energy along x at those temperatures. For weak 
disorder V — It as in (a), A T = A xx (q y — > 0) < K x indicat- 
ing the development of a finite superfluid stiffness D a from 
Eq. [m] with decreasing T. For strong disorder V — it as in 
(b), D s = at all T. 

In a system with a broken gauge symmetry, the longi- 
tudinal and transverse responses are no longer equal and 
their difference is precisely the superfluid stiffness D s or 
the related quantity, superfluid density p Sl given by 



p s = D s /tt 



[A L 
[K x 



A T ] 
-A T 1 



(14) 



It can be seen from Eq. |l4| that on a lattice the superfluid 
density at T — is indeed bounded above by the kinetic 
energy. In recent work |13] we have obtained an improved 



upper bound on D s in a disordered system in terms of 
the local kinetic energy which highlights the dominance 
of the weak links in determining the superfluid stiffness. 

In order to extract the superfluid stiffness D s from 
Eq. [h| we must extrapolate A xx (q y ) to q y — * 0. Using 
general symmetry arguments we have 



j Q (q) = A aj a(q)A J a(q) 

9a<7/3 \ . 1 



A Q /3 = I 5 a /3 



A-(^) + ^A-(^) (15) 



so that the linear term in the expansion of A T and A L 
is absent and the lowest order term is quadratic in q y . 
However, the momentum discretization on an 8 x 8 lattice 
is too coarse to see this quadratic behavior. 
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FIG. 8. The superfluid stiffness D s and K x as a function 
of T for V = 1.0, U = -it and (n) = 0.875. Also shown is 
the charge stiffness D at the lowest T = 0.1. 

It is clear from Fig. ^ that the transverse correlations 
behave quite differently from the longitudinal correla- 
tions. For weak disorder, at high temperature, A T ap- 
proaches K x , but as T is decreased, the two quantities 
no longer match, indicating that a nonzero superfluid 
density is developing as shown in Fig. [| We sec that 
D s becomes significantly different from zero at temper- 
atures T < 0.2i. This is consistent with estimates pq| 
which put T c m O.lt based on a finite size scaling anal- 
ysis of the pairing correlations, but seems to contradict 
recent suggestions that T c is much lower, approximately 
0.03i. In Fig. || we also show the behavior of K x which 
shows no special features as T is lowered. K x declines 
from 0.68 at V — to 0.39t at V = 5t, while D s changes 
by almost two orders of magnitude. While a reduction 
in hopping is expected in the presence of disorder, the 
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smooth behavior of the kinetic energy emphasizes that 
such a local quantity cannot serve as an order parameter 
for the localization transition. When disorder is strong, 
remains pinned at K x , for all T, suggesting that a 
superconducting phase is not present. 

Thus from the raw data itself there is compelling ev- 
idence for a superconducting phase at low temperature 
and at low disorder that is qualitatively distinct from the 
non- superconducting phase at higher disorder. 

Finally, we note that the mean field gap is of the or- 
der of the hopping integral t for U = —At, therefore 
quasiparticle excitations across the gap are suppressed 
by a factor ~ exp(— t/T) = cxp(— 10) at a temperature 
T = Q.lt. The finite temperature transition is thus dom- 
inated largely by thermal phase fluctuations. 



C. Superconductor-Insulator Transition 

In order to determine the location of the transition, 
we now present data at a set of disorder values which 
sweeps through the values V = 1 - 4 which we argued in 
the preceding section brackets the transition. In Fig. ^ 
we show the extrapolated values of A xx (q y ) and K x as 
a function of disorder. It is evident that the transition 
is driven by the variation of A T . In Fig. [To] we show D s 
as a function of disorder strength at fixed temperature 
T = O.lt, for U = -3t and U = -At. The decrease in D s 
with increasing disorder is consistent with the decline in 
the order parameter shown in Fig. ||. 



00 




V 



FIG. 9. The transverse current-current correlation func- 
tion A T = \\vci qy ^(, K xx {q y ) as a function of disorder at 
T — 0.10. Also shown is K x = A L , the longitudinal response 
function, as a function of disorder. The difference between 
A L and A T is the superfluid stiffness as seen from Eq. 
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FIG. 10. The superfluid stiffness D s and the charge stiff- 
ness D as a function of disorder strength V for U = — 3t and 
U = —At. Note the rapid suppression with disorder and the 
transition from a superconductor to an insulator beyond a 
critical disorder. 

The superfluid stiffness D s ~ 6^ where 5 = \V — 
Vc\/\Vc\ is the distance from critical disorder. The ex- 
ponent £ is expected to be larger than unity since £ = zv 
and in 2d it has been argued that v > 2/d = 1 and z = 2. 
A value of £ > 1 implies that the finite size rounding will 
shift the critical point on the infinite lattice to higher val- 
ues compared to the point where D s becomes small on 
finite lattices. So we expect that the critical point for the 
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SIT may lie around V c « 3-4t for U = -At. 
6 



Ld 



O 



O 




in the pair correlations and the transition to insulating 
behavior in the resistivity might reflect changes in the 
noninteracting eigenstates of the Hamiltonian. Is the fact 
that the pairing correlations are robust at V = but zero 
at V — 5t a consequence of some changes in the extent 
of the single-particle wavefunctions due to disorder? 

In Fig. [ll] we show the density of states N(E) for U = 
and different amounts of randomness bracketing V c . We 
see that disorder broadens N(E), as expected, but the 
behavior of this quantity through V c is smooth. 

We show in Fig. [l2| the localization length or the "size" 
of the eigenstate at the Fermi surface, defined by ^i oc — 
as a function of disorder strength. £j oc shows 
a smooth decrease as a function of V, without any sharp 
feature at V c . We conclude that the SIT is not occuring 
as a consequence of a U = Anderson transition on the 
finite lattice, even though the wavefunctions are localized 
on the scale of the linear lattice size L. Instead, the 
transition is a genuinely nontrivial many body effect. 



D. Coherence length 



FIG. 11. The density of states N(E) as a function of 
energy E for noninteracting fermions (U — 0) on an 8 x 8 
lattice at density (n) — 0.875 for disorder strengths around 
V = 2t, it < V c and V = At > V c , where V c is the critical 
disorder of the interacting problem with U = —At. 
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FIG. 12. The approximate localization length of the eigen- 
state at the Fermi surface inferred from the participation ratio 
by as a function of disorder strength V. We see that the single 
particle eigenstates do not show any sharp behavior around 
the critical disorder V c ~ 3.25 found for the SIT in our QMC 
simulations of the interacting problem. 

It is reasonable to ask to what extent the sharp drop 



In principle, we can extract the superconducting co- 
herence length £ for the many body problem from the 
dependence of A T (q y ) = a + bq 2 y for small q y . From Eq. |lj 
we see that 



D s (q y ) _ D s 



(16) 



where £ 2 = bj (D s /tt). As a function of disorder £ is found 
to decrease slightly for low disorder and is expected to 
diverge as the critical disorder is approached. However, 
it is difficult to deduce such a divergence from the data 
since both b and D s are becoming small near the transi- 
tion. Further work on this problem is required, since it 
would be useful to obtain the coherence length to track 
the quantum phase transition. 



VII. CHARGE STIFFNESS 

A superconductor is characterized by the Meissner ef- 
fect, measured by the superfluid stiffness, as well as by 
an infinite conductivity. A signature of the latter is a 
delta function in the optical conductivity 



Rea(oj) = D5(oj) +Rea re g(uj) 



(17) 



with weight D = n[K x — lim^^o R- e A X2 ,(q = 0;u> + i0 + )], 
known as the charge stiffness. The regular part of the 
conductivity is given by (suppressing the q = and omit- 
ting the xx subscripts) 



Recr reg (w) = 



ImA(w) 



(18) 



where A(w + i0 + ) = ReA(cj) + ilmA(w). 
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In order to obtain the dc limit we proceed as follows. 
We start with the sum rule 



f°° 7T 

/ dujRea(uj) = —K x 
Jo 2 

and combine with Eq. ^ to get 

J dujRea reg (uj) = -K x - —. 

Next, using the spectral representation for A(z), 

r duj ImA(w) 
A(«) = / — 

and substituting z = iw„ we get, 

... 2 /•°° J wImA(w) 

Using Eq. 



A(w„) = 



7T 







dwRecr re g(w) 



Reer re g(w) 



2 2 r / 

-w„ / 

7T ./ W Z + UT 



(19) 

(20) 

(21) 

(22) 

(23) 
(24) 



ii 



Substituting for the first term from Eq. 20 and defining 
the Matsubara correlation function 



D(w n ) = ir[K x - A(w„)] 



whence, 



poo 

Jo 



(25) 



(26) 



The behavior of A(w„) as a function of n is shown 
in Fig. |l3| for low disorder V = It in (a) and for high 
disorder V = At in (b). The behavior of A(ui n ) is qual- 
itatively similar to Fig. [?]. That is, at strong disorder, 
A(to n — > 0) ph K x at all temperatures and according to 
Eq. 25 this implies the charge stiffness D m 0, as is the 
superfluid stiffness D s . At weak disorder and at low T 
on the other hand, Muj n — > 0) < K x implying that D 
is nonzero. In Fig. 04 we show D(cu n ) as a function of 
n which is found to increase monotonically with n from 
D(uj n — ¥ 0) = D to D(u! n — > oo) = 7r(— i^a;) (not shown 
in the figure) but verified in the data. 

The behavior of D as a function of disorder extracted 
from Eq. ^ is shown in Fig. |l^. We see that D and D 8 
are within 10-20% of each other for all the parameters 
shown. Thus there is remarkable consistency between 
the superfluid stiffness D s and the strength D of the delta 
function in the optical conductivity, obtained from two 
very different correlation functions. 

Do these techniques give sensible results in the non- 
interacting, clean limit? For U = V = we find the 



charge stiffness D/tt — 0.79 = (—K x ) whereas the super- 
fluid stiffness D s /n = 0.0243 for filling (n) = 0.86 and 
T = O.lt on an 8 x 8 system. Thus our numerics are 
correctly telling us that free fcrmions are metallic with a 
nonzero D, but a very small D s , which will go to zero as 
the system size increases. 
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FIG. 13. The behavior of A(cj„) as a function of n at 
temperatures T — 0.17 and 0.10. The corresponding filled 
points at n = are the values of K x at those temperatures. 
For weak disorder V = t lim^^o A(tj n ) < K x indicating a 
finite charge stiffness D whereas at strong disorder V = 4t, 
lim^„^o A(oj n ) m K x implying that D = 0. 

While the approximate equality of D and D s in Fig. [l(] 
for a superconductor is a good check on the calculation, it 
emphasizes that the charge stiffness D at T = cannot 
be used to characterize the non-superconducting state 
for V > V c since neither dirty metals nor insulators have 
a (5 function in <j(u>) at lu = 0. Hence we turn to the 
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conductivity. 
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FIG. 14. The behavior of Z>(w„) defined in Eq. (|2J) as 
a function of re for T = 0.1, U = -4t and (re) = 0.875 for 
disorder strengths V = 1, 2, 3.25, 4, 5t. The corresponding 
filled points at ui n — are the extrapolated values of the 
charge stiffness D at those V. The critical disorder V c = 3.25 
is identified by the vanishing of D{uj n ). The straight line is 
a linear fit to the low u> n data whose slope is proportional to 
the critical conductivity at the transition from Eqj3C| 



VIII. CONDUCTIVITY 



The dc conductivity <Jdc = lim^^oReoVeg^) defined 
in Eq. [l8] is of considerable theoretical and experimental 
interest as it distinguishes the two non-superconducting 
phases- metal (above T c ) vs insulator. The fluctuation- 
dissipation theorem relates ImA(tj) which is required for 
the calculation of cr^c, to A(r) which is obtained from 
QMC data by 



A(r) 



du> exp(— lot) 
7r [1 — exp(— (3u))] 



ImA(w) 



(27) 



valid for < t < j3. However, the evaluation of ImA(w) 
requires an analytic continuation of noisy imaginary time 
data [?] which is difficult. We derive below an approxi- 
mate expression for a^ c p2fl , analogous to that introduced 
previously for the susceptibility P4fl , by noting that if one 
sets t = /3/2, the kernel in Eq. E7fcuts off contributions 
from high frequencies, and the important range of lu is 
restricted to increasingly small values as f3 becomes large. 
Therefore, at low enough temperatures one might replace 



O 

a 

\ LO h 

a 



o 



o 

CN 



o 




Q I I I I I I I I I I I I I I I l_ 

0.2 0.4 T 0.6 0.8 



o 

CM 



LO 



a 
a 

a o 



o 



o 
to 




o 

0.2 0.4 0.6 0.8 
T 

FIG. 15. The behavior of the resistivity p obtained from 
Eq. |28| as a function of T for various disorder strengths. Figs, 
correspond to U — —3, —4, —6 respectively. 



11 



ImA(u;) ~ wCTdc over the entire range of integration, 
which leads to the result 



O 



Ode 



/3 2 A(r = (3/2) 



(28) 



Note that Eq. £8| is only valid in the normal state (T c as 
O.lt) where lmA(uj) ~ wa& c at low frequencies. We will 
present a number of self-consistent checks of Eq. 



28 



the metallic state above T c of the superconductor and the 
localized phase. We defer a discussion of the extraction 
of the conductivity at the SIT to the next section. 

If Fig. [l5| we show the behavior of the resistivity p — 
1 /<7dc obtained from Eq. ^8] as a function of temperature. 
The resistivity shows a behavior qualitatively similar to 
that seen in experiment: when the control parameter, in 
this case disorder, is weak, the behavior is metallic and p 
decreases as T decreases. On the other hand, for strong 
disorder, the behavior is insulating and p increases as T 
decreases. Our plots are qualitatively similar to those 
observed experimentally, though the experimental range 
of resistivities is much greater. 

As is often done experimentally, data for p{T) at dif- 
ferent V can be replotted to show p(V) for different tem- 
peratures. For V < V c the resistivity decreases as T is 
lowered, while for V > V c the resistivity increases as T is 
lowered. This leads to a characteristic crossing pattern 
in p(V) which allows for an estimation of the critical 
amount of disorder V c as well as the critical resistance 
p{V c ) at the transition. Note that the crossing pattern 
does not follow from any deep scaling principle. Instead, 
it is merely a consequence of the monotonicity of the plots 
of p(T) for a given V, which, to within error bars, either 
steadily increase or decrease as T is changed. 

From Fig. || and Fig. we see clear evidence for a 
SIT at a critical disorder V C (U) whose dependence on the 
strength of the attraction is shown in Fig. n~7l 



It has recently been emphasized by Sachdev |45| that 
using Eq. ^8] to extract the resistivity is not applicable 
near a quantum phase transition as there is no scale in 
the problem. Note that it was assumed in the derivation 
of Eq. that below some scale which was independent of 
T, it was possible to assume that ImA(w) ~ tuo~dc- This 
assumption breaks down near a quantum critical point 
since by definition all scales become soft. Away from 
the transition, Eq. ^ gives a good description of Pdc(T), 
however, close to the transition, it cannot be used to 
extract the critical conductivity. The agreement of the 
transition point obtained by the conductivity crossing 
plots and the measurements of the superfluid and charge 
stiffness suggest that Eq. |28| has a useful range of validity. 

We discuss another potential method to extract the 
conductivity at the critical point. As seen in Fig. 14 at 
a critical disorder D vanishes. At this disorder assume 
that Recr(o;) — > Co = const, for frequencies u> < u> c , a 
cut-off value. 
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FIG. 16. The behavior of the resistivity p obtained from 
Eq. ^ as a function of V for various temperatures. Figs, 
correspond to U = —3, —4, —6 respectively. 
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Then from Eq. BQ 



D(w n ) = 2uj\ 



dw- 



5 M 



2cr |w n | tan 1 ( — 



2ui 



cL;^H (29) 



which in the limit of small Matsubara frequencies is given 
by 



D(iO n ) = 7r<7 |u; n | + O(0J n ) 



(30) 
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FIG. 17. The locus of critical disorder V C (U) for intermedi- 
ate couplings in the disorder V - attraction \U\ plane. The V c 
values are obtained by two independent methods, the vanish- 
ing of superfluid density D s (open circles) , and the crossing of 
the resistivity p (open squares) . The filled circle at U = V = 
emphasizes that all noninteracting states are localized for any 
non-zero disorder V in two dimensions. 

The conductivity at the critical point obtained from 
Eq. [SO] and from the crossing of the resistivity curves 
described above are in agreement to about 10 %. 
Near a quantum critical point, we expect 



a reg (uj,T, V = V c ) = a Q a(uj/T) 
From Eq. |2^, this implies that 
D(u n ) D(T) 



(31) 



T 



T 

G(T)+F(n) 



pOO 

87r 2 n 2 (jQ / dx- 

J x c 



(32) 



where x — lj/T. Thus D(ui n )/T is a sum of two terms- 
the first one G{T) = D(T) /T is only a function of T and 
the second term F(n) is only a function of n, with F(n — + 
0) = 0. We set V = V c ~ 3.25i and by extrapolating 
the behavior of D{w n )/T to n -> obtain G(T). In 
Fig. [l8|, we show the behavior of F(n) vs n — u) n /2nT at 



the critical point for various temperatures. The data are 
not found to scale, unlike our expectations at a critical 
point. Instead if we plot D{uj n ) vs uj n we see a remarkable 
scaling behavior of the data for various temperatures as 
seen in Fig. |l^. It is not really clear as to why the data 
when plotted as in Fig. Il8] does not scale. 
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18. The behavior of F(n) = D(lo„,T)/T 

0,T)/T as a function of n = u> n /(2-KT) defined 
at the critical disorder V c ~ 3.25f and various T. 
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FIG. 19. The behavior of D(ui„,T) as a function of uj n at 
the critical disorder V c ~ 3.25t and various T. 

It has been claimed in Ref. [Q that since in QMC the 
lowest frequency that can be accessed is lu± — 2irT > T, 
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it is not possible to extract the dc resistivity using Eq. |30j 
in the low frequency limit. While this objection appears 
very sound, it is nevertheless the case that the conductiv- 
ity inferred from Eq. ^o| including its values in the vicinity 
of the critical point, is consistent with many other, com- 
pletely rigorously founded, aspects of our simulation. By 
this we mean that the location of the transition inferred 
from the analysis of the data using Eq. |3^ is in remarkable 
agreement with the location obtained from the superfluid 
stiffness D s , and the charge stiffness D. Furthermore, 
the value of the conductivity at the transition is consis- 
tent with the value obtained from Eq. ^8|. At present we 
do not understand fully why the method appears to be 
so consistent with our other data despite the objections 
raised in Ref. @. 

IX. CONCLUSIONS 

We have studied the effect of disorder on an s-wave su- 
perconductor of fixed coupling strength (modeled as an 
attractive Hubbard model away from half filling). We 
have found that with increasing disorder, the superfluid 
stiffness (obtained from the transverse current-current 
correlation function) and the charge stiffness (obtained 
from the r— dependent current-current correlation func- 
tion, vanish at a critical disorder, signaling a transition 
to a localized phase. The importance of our work lies in 
the fact that the SIT which has been observed experimen- 
tally, has eluded all mean field treatments of the problem. 
Ours is the first theoretical study of a fermionic model 
to obtain a transition between the superconductor and 
localized phases upon increasing the disorder strength. 



Once the location and exponents characterizing the 
transition are determined, the key question is the value of 
the resistivity at the transition and the possibility of its 
universality. There is some experimental evidence that 
despite the wide range of materials and control param- 
eters, the value of the resistance right at the transition 
R* is always quite close to the "universal" value |^,||,[5],|} 
Rq = h/Ae 2 . While there is still some debate concerning 
whether this number is truly the same for all systems, it is 
certainly clear that the variation in R* is much less than 
the variations in the location of the transition in other 
control parameters such as the temperature, magnetic 
field strength, or film thickness. Recent experiments of 
Yazdani et al. || have interpreted the variation in R* 
that exists in terms of separate bosonic and fermionic 
contributions to the resistivity. Thus, calculations with 
models which include electronic degrees of freedom like 
the attractive Hubbard Hamiltonian are needed to sup- 
plement work on bosonic theories. To address this set 
of issues concerning R* , we require an exact method to 
calculate the resistivity at the transition, as would be 
provided by maximum entropy techniques. We are cur- 
rently working on this problem. 
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X. OUTSTANDING QUESTIONS 

Having established the existence of the SIT the out- 
standing questions now relate to obtaining a quantitative 
characterization of the transition. For this it is necessary 
to perform finite size scaling in both the spatial (L — > oo) 
and the temporal (T — » 0) dimensions to obtain the lo- 
cation of the critical disorder from the vanishing of the 
superfluid stiffness D s as well as the vanishing of the 
charge stiffness D. From the scaling of the data it is 
then possible to extract the dynamical exponent z and 
the correlation length exponent v. Such an analysis will 
tell us whether the fermion SIT is in the same universality 
class as the bosonic superfluid-insulator transition or not. 
While there have been several studies of the SF-I transi- 
tion @H in the boson Hubbard model |g||(J[l8| and 
its variants Jl7j , we believe that the situation with regard 
to the value of the exponents is still unclear This 
is largely because of the complications of finite size scal- 
ing analysis inherent in a quantum phase transition that 
necessarily involves two variables-(system size L — > oo 
and temperature T — > Q). 
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